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ARSTRACT: 

The problem of preconditioning the matrices arising from pseudo-spectral Chebyshev 
approximations of first order operators is considered in both one and two dimensions. In 
one dimension a preconditioner represented by a full matirx which leads to preconditioned 
eigenvalues that are real, positive and lie between 1 and ill 2, is already available. Since 
there are cases in which it is not computationally convenient to work with such a 
preconditioner, we study a large number of preconditioners which are "more sparse" (in 
particular three and four diagonal matrices). The eigenvalues of such preconditioned 
matrices are compared. In particular, the analysis is carried out for the quantity 
max I Xi I /min I Xi I , where X ; are the preconditioned eigenvalues. 

We apply the results to the problem of finding the steady state solution to an 
equation of the type u t = + f, where Chebyshev collocation is used for the spatial 

variable and time discretization is performed by the Richardson method. 

In two dimensions different preconditioners are proposed for the matrix which arises 
from the pseudo-spectral discretization of the steady state problem in the square 
A = {(x,y,): - l«x<U, -KySl} 

U x + U y = f 
U(x, y, 0) = U 0 

with boundary conditions at x = 1 and y = 1. Results are given for the CPU time and 
the number of iterations using a Richardson iteration method for the unpreconditioned 
and preconditioned cases. 
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1. INTRODUCTION 

To obtain the pseudo-spectral or collocation approximation let P N be an interpolation 
operator. Let f(x) be a sufficiently smooth function defined in [-1,1] where f(x)=0 at the 
appropriate boundaries which yields a well-posed problem for (1.1). Then P N f is the 
interpolation of f at the collocation points Xj, i.e. 

P N f(Xj) = f(Xj) and P N f€B N . j = 0 N 

To obtain a Chebyshev Gauss-Lobatto pseudo-spectal approximation in the interval 
[-1,1] we choose Xj = cos jn/N (j = 0,...,N), which when j t 0,N are the extrema of the 
N th order Chebyshev polynomials T N (x) = cos(N cos _1 x). In order to construct the 
interpolant of f(x) at x, we define the polynomials 

(T x 2 )T^j(x)(-iy + 1 

gj(x) « G = 0 N) 

Cj N 2 (x - x^ 

C 0 = c^ = 2, Cj = 1 (lSjSN-1). 

One can easily see that gj(x k ) = 6j k . 

The N th degree interpolation polynomial P N f to f is given by 
N 

(1.8) P N f(x) = I f(Xj)gj(x) xeIR 

j=o 

We must now be able to express derivatives of P N f in terms of f at the collocation 
points Xj. Differentiating (1.8) we obtain 

d n P N f(x) N d n 

(1.9) = l f()d) g (x) 

dx 11 j=o dx n 

so that 

d n P N f(x k ) N 

= I f(xp (D n ) kj 

dx n j=o 


( 1 . 10 ) 
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where 

(' ll) (°nV ■ S~ 8 k w 

X = Xj 

The pseudo-spectral Chebyshev derivative operator can be represented by the N x N 
matrix S„ = [s^, 
where 

Cj(-l)j +k 

s jk = 

c k (Xj - X k ) 

• x j 

s jj = , 

2(1 - Xj 2 ) 


(k * )) 

2N 2 + 1 

S 00 " g = * S NN 


In particular the Chebyshev pseudo-spectral approximation for u = u , u(x,0) = u o (x) 
N 

is given by u N = I u(Xj,t)gj(x) and 


k=0 


9u 


N 


N 


(Xj.t) = I u(Xj,t)s k j . 
k=0 


at 


F 

t 
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2. OT ITLINE OF THE PROBLEM 
In the interval [-1,1] let 


2 j - 1 

L = cos — — n O' = 1.2,..., N) 
1 2N 

lies between Xj and Xj.j 



The modes Xj and ^ have the following properties: 

T N «p = 0 j = UN, and T N (x p = (-1) 1 i = 0 N 

Then consider the pseudo-spectral Chebyshev derivative operator with homogeneous 
boundary conditions at x = 1. This operator can be represented by the N * N matrix 



if i = j = N 


if i * j 


if i = j = 1,...,N-1 


The matrix S N is full. The condition number C(S N ) of S N is large. We have the 

following result which was obtained by Daniele Funaro. 

Lemma 1-1 : The condition number of S N increases at least like N 2 . 

Proof : Let 0 • U denote the norm in £(IR^JR n ). Then the condition number 

of S N is given by 

c(S N ) = IS N H IS^B . 


( 2 . 1 ) 
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It is known that I S N I >p(S N ) ^CjN 2 , where p(S N ) is the spectral radius of S N and 
Cj is a constant independent of N. On the other hand, we have 


( 2 . 2 ) 


IISJ/8 = sup !S- N V R N. 
I?l = 1 

R n 


Choose 


% = 


i/n 

j = Furthermore, 


Then II <Po H ^ - 1 and (S^<po)j - — ( x j 


(2-3) 


4-KV !) 2 > n ~ { ( x ' !) 2 wdx = 

j- 1 -l 


1), for 


where c 2 does not depend on N. 

This implies 8 S^j 1 It > B Sg 1 ^ ^ n >c 2 - Finally, using (2.1), we get C(S n )^C 3 N 2 . This 

proves the claim. 

Although the condition number is particularly meaningful for numerical applications, 

its determination is generally very difficult. Another quantity which is meaningful for 

practical application is o(M) = max I 1 /min 1 X- 1 where M is an N x N matrix and Xj 

i i 1 

i = 1,...,N are its eigenvalues. It can be shown empirically that a(S N ) behaves like N. 

We are interested in finding a "preconditioner" for S N . In particular we are 
concerned with finding a matrix R N such that the quantitity a(R^ 1 S N ) is small. In 
general this does not imply that the corresponding condition number will be small (so the 
word "preconditioner" is not correct). 

In [DF ] S N is preconditioned by R N = where the N x N matrix D N = {d^} 

is defined by 

d ii = '^i-T x i> > = 1 N 

d ii-l = l/(x,_ r Xj) i = 2....N 
djj = 0 otherwise 
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Hence D is the upwind finite differences matrix relative to the grid x j. is 

the operator which maps the values of a polynomial in P N-1 at the staggered grid points 
{$!>-, 5 N } into the values at the mesh points {x 1 ,...^c N }. 

Preconditioning by Z N D N results in preconditioned eigenvalues that are real, positive 
and lie between 1 and 77/2. The ratio o(D^ f 1 Zj v , 1 S N ) is bounded by fit 2 (see [DF]). 
This is particularly interesting when the solution of the system 
(2.4) S n u + f = 0 

has to be found. If an iterative method is used, iterating M N = (Z N D N ) _1 S N instead of 
S N results in convergence in a few iterations. In the end of the computation the system 

(Z n D n )u + f = 0 has to be solved. So we require that the matrix R N = Z N D N can be 

inverted easily. Although Z is a full matrix, it can be inverted very inexpensively in N 

log N operations. Thus the matrix R N = Z^D^ can be inverted very inexpensively. 

Nevertheless there are cases in which it is not computationally convenient to work with a 
full preconditioner. 

In particular when using an implicit method to find the solution at time T > 0 of 
the equation 


u t = S N u + f, 

(2.5) u(x,0) = u Q (x), 

u(l,t) = 0 0<StST. 

If the implicit Euler method is used, iterates of the matrix (I + AtS N ) are considered. A 
good preconditioner for this matrix turns out to be the matrix (I + AtZ^D^,). 
Unfortunately, due to the fact that Zj^ is a full matrix, (I + AtZ N D N ) cannot be inverted 
inexpensively. In this work we shall present a large number of preconditioners which can 
be applied to the situations illustrated above. 
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.3. ANALYSIS IN ONE DIMENSION 

In order to have the matrix I +At Z N D N which can be easily inverted we substitute 
Z N by some suitable matrix which has a simpler form and which has to be regarded as 
an approximation of the operator related to Z N . 

The first idea is the following. Take Z N = {zy} such that 


A/ 

Z ii 

= .5 

Z ii + 1 

= .5 

A/ 


Z NN 

= 


Then Z N D N is a tridiagonal matrix so tha I + At Z N D N is also a tridiagonal matrix. It 

A/ 

can be shown that the eigenvalues of M N = (Z N D N ) -1 S N take the following form. 


k-l 

\ = k I Xj 


k = 


Thus, the eigenvalues of the preconditioned matrix are real and positive. Hence we 
have o(M n ) = N. The choice of as in (3.1) corresponds to shifting the values from 
the staggered mesh to the initial mesh by averaging the two neighbour values. Instead of 


A/ A/ 

this we can choose Z N = {Zy} corresponding to interpolation by first order polynomials. 

A* A/ 


Tliis leads to the following definition 
x i'S i+ i 


z ii = 


Si - Si 


i+l 


of the matrix Z N = {Zy} 
i = 1,...,N-1 


( 3 . 2 ) Xj - (j 

z ihl = -T T * = UJ-1 

i ’i+l m 

We found emperically, that the eigenvalues of the preconditioned matrix are still real and 
positive, but the quantity o(M N ) is now worse than that of the previous case. 

Another simple preconditioner which this time is not of the form Z N D N is defined 
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by R n = { rjj } where 

Tjj = 0 i = 

r NN = ^( X N-1 ‘ %) 

(3.3) r i i+1 = 1%! - x j+1 ) i = 1.....N-1 

r i,i-i = ^-l ‘ x i + i> > = 2,...,N-1 

r N-l,N = ‘ X n) 


In this case we still get real and positive preconditioned eigenvalues. Namely, they take 
the form: 


X„ = 


k sin tt/N 


k = 1,...,N-1 


u{x N }, 


L ''i - x| J 

where X N is approximately equal to 2.46. We have o(M N ) = N-l. This is the best we 
tested using tridiagonal preconditioners. Up to now the improvements are poor so we 
have to consider better approximations of the matrix Z N . 

We consider an approximation of the operator related to by interpolation with a 
polynomial of degree 2. One possible choice is the following: Z = (z-} where 


N 

i— ‘ 

II 

(x x = 

= 5 2 )/«i 

' *>?) 





N > 
Si 
II 

( X 1 * 

■ ?i)/(5 2 

- U 





z U-i 

= (Xj 

- ti)(Xj - 

^iV((5i.i 

- - 


i = 2,.. 

.,N-1 

z ii = 

- 




i + i» 

i = 2,.. 

„N-1 

Z i,i+1 

= ( x i 

- ?i-i)( x i 


- ^i-lX^i+1 

- ^i)) 

i = 2,. 

,.,N-1 

Z N,N-1 

= (x 

N ‘ W^N-l ' 





A 

Z NN : 

= ( X N 

- Sn-I^N ' ?N-j) 






Now Z is a three-diagonal matrix, hence ZD is a four-diagonal matrix. The numerical 
experiments performed up to N = 32 give the following results. The eigenvalues X of 
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M n = (Z n D n )' 1 S n are in general complex. The have positive real part greater than .98. 
Most of them are concentrated at 1. Figure 3.1 shows the behavior of the X’s for N = 
12 and N = 20. 




Figure 3.1 - Location of the preconditioned eigenvalues in the complex plane. 

The corresponding o(M N ) is represented in Table 3.1 for various N. This time 
c(M n ) is bounded with respect to N. 


N 

°(M n ) 

8 

2.602 

16 

2.724 

24 

2.758 

32 ! 

2.770 


Table 3.1 - Case of Four-diagonal preconditioner 
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Other possibilities for the tridiagonal matrix Z were tested. For example one 

A 

possible choice, analogous to that of the matrix in (3.1), is to take Z N = {Zy} such that 

Vi - -V* 

(3-5) z H = 3/4 

= 3/8 

Among all the experiments, the matrix proposed in (3.4) gives the best results. 
Five-diagonal preconditioners were tested, through interpolation by third degree polynomials. 
The results do not imporve those corresponding to Table 3.1. 

Now, if the system of linear equations (1.5) is solved, for example, by the implicit 
Richardson method, we are concerned with o(M N ), where M N = (I + At Z N D N ) _1 (I + At S N ,) 
and Z N is the matrix in (3.4). The graph of o(M N ) versus At is reported in Figure 2.2 
for some values of N. So, in addition to the fast convergence of the iterative scheme, 
the preconditioning matrix can be inverted efficiently. 


* 



1 


2 3 4 


5 


_Afc 


Figure 3.2 



- 11 - 


Now, if the steady state solution of problem (2.5) has to be found (which 
corresponds to the solution of problem (2.4), the explicit Richardson method can be used. 
This leads to the iterative scheme 

(3.6) u n+1 = (I + AtS N )u n , n€ N 

If X denotes the general eigenvalue of S N , the scheme is stable provided we choose At 
such that 

f -2Re\ 1 

(3.7) 0 < At < inf . 

X 1 X | 2 J 

Within the stability region we have p(I + AtS N ) < 1 where p denotes the spectral radius. 

In order to speed up the convergence we experimentally find At* such that 
p(I + At*S N ) attains its minimum inside the interval of stability. In general At* is not 
available. We use it here only to compare preconditioners. The same experiments are 
made for the preconditioned schemes, i.e.: 

(3.8) u n+1 = (I + AtZ N D N )- 1 (I + AtS N )u n 

(3.9) u n+1 = (I + AtZ N D N )' 1 (I + AtS N )u n 

where Z N and Z N are respectively, the full matrix proposed in [DF] and the tridiagonal 
matrix given by (3.4). 

Both the schemes (3.8) and (3.9) converge to the same solution of (3.6). The 
results of these experiments are reported in Tables 3-2, 3-3, and 3-4 


N 

Maximum At for Stability 

At* 

p Corresponding to At* 

8 

.05461 

.02724 

.9817 

16 

.01836 

.009120 

.9919 

24 

.009232 

.005244 

.9939 

32 

.005138 

.002521 

.9966 


Table 3-2 Minimum spectral radius of the unpreconditioned amplification matrix 
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N 

Maximum At for Stability 

At* 

p Corresponding to At* 

8 

1.281 


.7810 

.2190 

16 

1.275 


.7787 

.2213 

24 

1.274 


.7783 

.2217 

32 

1.274 


.7782 

.2218 

Table 3-3 Minimum spectral radius 

of the preconditioned matrix: case of Z N . 

N 

| Maximum At for Stability 

At* 

p Corresponding to At* 

8 

.7685 


.5552 

.4448 

16 

.6381 


.3190 

.7112 

24 

.4258 


.2129 

.8465 

32 

.3210 


.1605 

.9043 


Table 3-4 Minimum spectral spectral radius of the preconditioned matrix : case of Z N r 
We also considered the second order Runge-Kutta scheme. This leads to the iterative 
scheme 

At^ 

(3.10) u n+1 = (I + AtS N + — — S 2 )u n , m-N 
In this case the stability restriction on At is given by 

(3.11) At 3 1 X j 4 + 4 At 2 (ReX) | X | 2 + 8At(ReX) 2 + 8ReX < 0. 


for all eigenvalues, X, of S N - We obtained results that were qualitatively analogous to 
that of Tables 3-2, 3-3, and 3-4. 
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4. ANALYSIS IN TWO DIMENSIONS 

In two dimensions we will consider the following steady state problem on the square 
A = {(x,y): -H x ( 1, -1 t y U} with homogeneous boundary conditions at x = 1, 
and y = 1. 

(4 A) u x + Uy = f 

u(x,y,0) = u 0 (x,y)€A 

u(l,y,t) = u(x,l,t) = 0. t £ 0 

The essential idea in obtaining a pseudo-spectral approximation to (4.1) is the same as it 
was in 1-dimension. That is, to approximate spacial derivatives by constructing a global 
interpolant through discrete points. To obtain the Chebyshev pseudo-spectral approximation 
we take as these points Xj = y- = cos zzj/N for j = This means that we must 

interpolate at the N 2 points (Xj,yp for i,j = 1,2,...,N. Consequently, the Chebyshev 

derivative operator for this problem can be represented by the N 2 x N 2 matrix + 

P £ S n < 2 >P, where S^ 2 > is a block diagonal matrix whose blocks are each equal to S N , and P 
is a permutation matrix. If one orders the N 2 points (Xj,yp by rows then S^ 2 -* 
corresponds to the derivative in the x-direction and P is constructed so that P t S^, 2 ^P 
corresponds to the derivative in the y-direction. Without preconditioning + P c SL 2 ^P 
is ill-conditioned. 

As we saw in section 3 ZD is a good preconditioner for S N . Thus, a natural 

approach to finding a preconditioner for + P £ Sj^P is to try + P'Z^D^P, 

where and D® are N 2 x N 2 block diagonal matrices whose blocks are the N * N 
matrices Z and D, respectively. 

To analyze the behavior of the eigenvalues of the preconditioning matrix, we define V 
as the generic eigenvalue of the preconditioned matrix, p N = max 1 Xj | /min I V | (i =1,2,.. .N 2 ), 
o N is the maximum o such that ReX^o, and r N is the minimum r such that | X- 1 1 «r. 
In particular, r N and o N give us an idea of the location of the eigenvalues. (See figure 
4-1 below.) 
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Numerical experiments performed for N = 4,6,8 are summarized in Table 4-1 



Table 4-1 

Although the eigenvalues of this preconditioned matrix, (Z^D^ + P t Z^D^P)' l (S^ + 
P 1 S< n 2 >P), are well behaved the matrix is full and thus difficult to invert. Another 
approach to constructing a preconditioner is to substitute in Z^D^ + P l Z^D^P the 
tridiagonal matrix Z defined by (2.4) in place of Z. We will denote this new N 2 x N^ 2 -* 

A A 

matrix by Z^D^ + P l Z^D^P. This matrix represents a finite difference scheme 
depending on seven points as illustrated by the stencil in Figure 4-2. 

• 0 + lj) 

• • • « 

(i,j-2) (i,j-l) (i,j) (i,j + 1) 

• C-lj) 

• (i-2j) 


Figure 4-2 
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Results similar to those presented in Table 4-1 are prsented in Table 4-2 for 
z( 2 >d ( 2) + P l Z^ 2) D (2) P. 


N 

p N 

r N 

°N 

4 

2.206 

.980 

.897 

6 

5.026 

1.453 

.488 

8 

8.494 

1.602 

.306 


Table 4-2 

Although p N corresponding to Z®D^ + P 1 Z^D^P increases more quickly than p N , 
corresponding to Z^D® + P l Z^D^P, the matrix can be inverted more efficiently. This 
is because Z^D^ + P l Z^D^P is a banded matrix (with N lower codiagonals and 2N 
upper codiagonals). 

Another preconditioner that we considered was of the form 

+ P t D (2) P). 

As in the 1-dimensional case, if the steady-state solution is to be found, the explicit 
Richardson method can be used. We experimentally find At* such that p(I + At*W N ) 
attains its minimum inside the region of stability. The same experiments are made for the 
preconditioned matrices. The results of the experiments are reported in Tables 3-3, 3-4, 
and 3-5. 


N 

Maximum At for Stability 

At* 

p Corresponding to At* 

4 

.1509563 

.07547812 

.9248009 

6 

.05228855 

.02614427 

.9715761 

8 

.0273053 

.01365251 

.9817108 

10 

.01944351 

.009721750 

.9810511 


Table 4-3 Minimum spectral radius of the unpreconditioned amplification matrix 
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N 

Maximum At for Stability 

At* 

p Corresponding to At* 

4 

.9879928 

.7172827 

.6693591 

6 

.9094622 

.6584506 

.8027519 

8 

.7200746 

.4954113 

.8715243 

0 

.6038381 

.3623029 

.8995660 


Table 4-4 Minimum spectral radius of the preconditioned matrix: 


Case of Z N P l Z n P(D n + P l D N P) 


N 

Maxiumum At for Stability 

At* 

p Corresponding to At* 

4 

1.009613 

.6946138 

.3763489 

6 

.815170 

.6798409 

.6681257 

8 

.7685055 

.6870439 

37895085 


Table 4-5 Minimum spectral radius of the preconditioned matrix: 

Case of ZjPn + P l Z N ,D N P 

We also considered the second order Runge-Kutta scheme, and we obtained results 
that were qualitatively analogous to those of Tables 4-3, 4-4, and 4-5. 

We applied the Richardson schemes in the unpreconditioned version and in the 

A A /V 

preconditioned versions using the preconditioners Z N D N + P 1 Z N D N P and 

Z N P t Z N P(D N + P^^jP), to find the solution of the model problem 

u t = + Uy - cxsin(a(x+l)) + asin(a(y+l) 

u(x,y,0) = sin(y-l)sin(x-l) 
u(-l,y,t) = u(x,-l,t) = 0 

of the type in (4.1). We used the optimal At* listed in tables 4-3, 4-4, and 4-5. Using 

the norm we considered the scheme to converge when the exact error stabilized. We 

also calculated global CPU times. The experiments were performed on the IBM 3081 and 
the results are reported in Tables 4-6, 4-7, and 4-8. 
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N 

No. of Iteration 
for Convergence 

Error 

CPU Time for Convergence 
(in seconds) 

4 

122 

8 

X 

H-* 

o 

1 

.21 

6 

400 

.104182 x 10' 3 

1.33 

8 

900 

.47112 x 10 -6 

5.82 


Table 4-6 Unpreconditional Euler 


N 

No. of Interation 
for Convergence 

Error 

CPU Time for 
Convergence 

CPU Time for Construction 
of Preconditioner and 
Finding LU Decomposition 
of Preconditioner 

4 

11 

.1246626 x 10 - 1 

.03 

.00 

6 

34 

.10400 x 10 ‘ 3 

.26 

.04 

8 

69 | 

.4709596 x 10' 6 

1.17 

.11 


Table 4-7 Preconditioned Richardson Case of Z^D + P l Z^DP 


N 

No. of Interations 
for Convergence 

Error 

CPU Time for 
Convergence 

CPU Time for 
Construction of 
Precondtioner and 
Finding LU 
Decomposition of 
Preconditioner 

4 

32 

.12468 x 10' 1 

.09 

.00 

6 

173 

.104255 x 10' 3 

1.10 

.00 

8 

263 

.4709613 x IQ- 6 

3.31 

.00 


Table 4-8 Precondition Richardson Case of Z N P l Z N P(D + P l DP). 
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